PHYSICAL  REVIEW  B  73,  184113  (2006) 


Time-resolved  diffraction  profiles  and  atomic  dynamics  in  short-pulse  laser-induced  structural 

transformations:  Molecular  dynamics  study 

Zhibin  Lin  and  Leonid  V.  Zhigilei* 

Department  of  Materials  Science  and  Engineering,  University  of  Virginia,  116  Engineer’s  Way,  Charlottesville,  Virginia  22904-4745,  USA 
(Received  22  October  2005;  revised  manuscript  received  20  January  2006;  published  16  May  2006) 


The  diffraction  profiles  and  density  correlation  functions  are  calculated  for  transient  atomic  configurations 
generated  in  molecular  dynamics  simulations  of  a  20  nm  Au  film  irradiated  with  200  fs  laser  pulses  of  different 
intensity.  The  results  of  the  calculations  provide  an  opportunity  to  directly  relate  the  detailed  information  on  the 
atomic-level  structural  rearrangements  available  from  the  simulations  to  the  diffraction  spectra  measured  in 
time-resolved  x-ray  and  electron  diffraction  experiments.  Three  processes  are  found  to  be  responsible  for  the 
evolution  of  the  diffraction  profiles.  During  the  first  several  picoseconds  after  the  laser  excitation,  the  decrease 
of  the  intensity  of  the  diffraction  peaks  is  largely  due  to  the  increasing  amplitude  of  thermal  atomic  vibrations 
and  can  be  well  described  by  the  Debye-Waller  factor.  The  effect  of  thermoelastic  deformation  of  the  film  prior 
to  melting  is  reflected  in  shifts  and  splittings  of  the  diffraction  peaks,  providing  an  opportunity  for  experimental 
probing  of  the  ultrafast  deformations.  Finally,  the  onset  of  the  melting  process  results  in  complete  disappear¬ 
ance  of  the  crystalline  diffraction  peaks.  The  homogeneous  nucleation  of  a  large  number  of  liquid  regions 
throughout  the  film  is  found  to  be  more  effective  in  reducing  long-range  correlations  in  atomic  positions  and 
diminishing  the  diffraction  peaks  as  compared  to  the  heterogeneous  melting  by  melting  front  propagation.  For 
the  same  fraction  of  atoms  retaining  the  local  crystalline  environment,  the  diffraction  peaks  are  more  pro¬ 
nounced  in  heterogeneous  melting.  A  detailed  analysis  of  the  real  space  correlations  in  atomic  positions  is  also 
performed  and  the  atomic-level  picture  behind  the  experimentally  observed  fast  disappearance  of  the  correla¬ 
tion  peak  corresponding  to  the  second  nearest  neighbors  in  the  fee  lattice  during  the  laser  heating  and  melting 
processes  is  revealed. 
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I.  INTRODUCTION 

Short  (pico-  and  femtosecond)  pulse  laser  irradiation  has 
the  ability  to  bring  material  into  a  highly  nonequilibrium 
state  and  provides  a  unique  opportunity  to  study  the  material 
behavior  and  phase  transition  dynamics  under  extreme  con¬ 
ditions.  The  challenge  of  probing  fast  structural  transforma¬ 
tions  is  being  met  by  active  development  of  a  variety  of  time 
resolved  probe  techniques. 1-1 6  Until  recently  most  of  the  data 
on  the  kinetics  of  laser-induced  phase  transformations  has 
been  provided  by  optical  probe  techniques,  e.g..  Refs.  1-4. 
While  high  temporal  resolution  is  readily  achievable  in  opti¬ 
cal  pump-probe  experiments,  the  reflectivity  measured  by 
optical  probes  can  only  reveal  changes  in  the  electronic 
structure  of  the  irradiated  surface  and  provides  limited  direct 
information  on  atomic  structural  rearrangements. 

Recent  advances  in  time-resolved  x-ray  and  electron  dif¬ 
fraction  techniques  open  up  an  exciting  opportunity  to  go 
beyond  the  analysis  of  the  characteristic  time  scales  of  laser- 
induced  phase  transformations  and  to  directly  probe  the  tran¬ 
sient  atomic  dynamics.  For  example,  observations  of  the  in¬ 
ertial  motion  of  atoms  on  the  optically  modified  or  softened 
potential  energy  landscape,  reported  by  Lindenberg  et  al.,6 
provide  new  insights  into  the  mechanisms  of  nonthermal 
melting  of  covalently  bonded  materials.  Excitation  of  large 
coherent  atomic  displacements  at  low  laser  fluences  and  dis¬ 
ordering  or  melting  at  higher  fluences  has  been  deduced  by 
Sokolowski-Tinten  et  al.  from  analysis  of  time  evolution  of 
the  diffraction  signals  obtained  for  50  nm  bismuth  films  ir¬ 
radiated  with  femtosecond  laser  pulses.7  The  resolidification 


process  of  a  laser-melted  surface  region  of  an  InSb  target  has 
been  studied  with  nanosecond  temporal  resolution  by  Harbst 
et  al.  and  the  velocity  of  the  resolidification  front  has  been 
measured  for  different  laser  fluences.8  Due  to  the  limited 
intensity  of  the  available  x-ray  pulses,  the  diffraction  signals 
are  typically  derived  from  rocking  curves,  for  fixed  Bragg 
angles.  The  advances  in  short-pulsed  electron  sources  enable 
a  competitive  alternative  to  x  rays  in  the  exploration  of 
atomic  dynamics.9  High  structural  sensitivity  and  subpico¬ 
second  time  resolution  were  recently  demonstrated  by  Si- 
wick  et  al.  in  an  electron  diffraction  study  of  ultrafast  solid- 
to-liquid  transition  dynamics  in  20  nm  aluminum  films 
irradiated  with  120  fs  laser  pulses.11  The  diffraction  intensity 
over  a  range  of  scattering  vectors  was  measured  in  this  work, 
allowing  for  analysis  of  time  evolution  of  the  density  corre¬ 
lation  function  during  the  melting  process  and  providing  in¬ 
formation  on  the  atomic  rearrangements  during  the  first  pi¬ 
coseconds  following  the  optical  excitation. 

Although  time-resolved  diffraction  experiments  provide 
important  atomic-level  insights  into  the  fast  laser-induced 
processes,  the  complexity  of  the  nonequilibrium  phase  trans¬ 
formations  hinders  the  direct  translation  of  the  diffraction 
profiles  to  the  transient  atomic  structures.  Atomic-level  simu¬ 
lations  can  help  in  reliable  interpretation  of  experimental 
observations.  Indeed,  classical  and  ab  initio  molecular  dy¬ 
namics  (MD)  simulations  have  been  used  to  study  the 
mechanisms  and  kinetics  of  laser-induced  nonthermal17'18 
and  thermal19-21  melting  processes,  the  evolution  of  voids  in 
photomechanical  spallation,22  as  well  as  the  dynamics  of  ex¬ 
plosive  material  disintegration  and  ablation.23-24  MD  simula- 
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tions  provide  complete  atomic-level  information  on  the 
mechanisms  of  laser-induced  phase  transformations.  At  the 
same  time,  diffraction  profiles  and  density  correlation  func¬ 
tions  can  be  calculated  from  atomic  configurations  predicted 
in  MD  simulations,  providing  a  direct  connection  between 
the  results  of  MD  simulations  and  time-resolved  diffraction 
experiments. 

In  this  paper  we  report  the  results  of  calculations  of  dif¬ 
fraction  profiles  and  density  correlation  functions  from  tran¬ 
sient  atomic  configurations  obtained  in  MD  simulations  of 
short-pulse  laser  melting  of  20  nm  Au  and  A1  films.  The 
computational  model  used  in  MD  simulations  of  laser  melt¬ 
ing  of  metal  films  is  briefly  described  next,  in  Sec.  II.  Nu¬ 
merical  methods  for  calculation  of  the  structure  functions 
from  atomic  configurations  are  discussed  in  Sec.  III.  The 
microscopic  picture  of  the  competition  between  the  homoge¬ 
neous  and  heterogeneous  melting  processes  obtained  in  MD 
simulations  is  discussed  and  related  to  the  evolution  of  the 
diffraction  profiles  and  density  correlation  functions  in  Sec. 
IV.  Connections  between  the  characteristic  features  of  the 
diffraction  profiles  and  the  mechanisms  of  laser  melting  re¬ 
vealed  in  the  simulations,  as  well  as  the  implications  of  the 
simulation  results  for  interpretation  of  experimental  data  are 
reviewed  in  Sec.  V. 

II.  COMPUTATIONAL  SETUP  FOR  SIMULATIONS  OF 
LASER  INTERACTIONS  WITH  Au  THIN  FILMS 

MD  simulations  of  fast  laser-induced  structural  transfor¬ 
mations  are  performed  for  thin,  20  nm,  freestanding  metal 
films  irradiated  by  a  short,  200  fs,  laser  pulse.  This  choice  of 
the  computational  system  is  defined  by  the  availability  of 
high-quality  time-resolved  electron  diffraction  data  obtained 
in  the  transmission  mode  for  thin  freestanding  films. 10-16 
Most  of  the  simulations  reported  in  this  paper  are  for  Au 
films,  with  an  additional  simulation  performed  for  an  A1  film 
with  irradiation  conditions  comparable  to  the  ones  used  in  a 
recent  experimental  study.11  The  simulations  are  performed 
with  a  hybrid  atomistic-continuum  model  that  combines  the 
classical  MD  method  for  simulation  of  nonequilibrium  pro¬ 
cesses  of  lattice  superheating,  deformation,  and  melting  with 
a  continuum  description  of  the  laser  excitation  and  subse¬ 
quent  relaxation  of  the  conduction  band  electrons.  The  model 
is  based  on  well-known  two-temperature  model25  (TTM), 
which  describes  the  time  evolution  of  the  lattice  and  electron 
temperatures  by  two  coupled  nonlinear  differential  equa¬ 
tions.  In  the  combined  TTM-MD  method,  MD  substitutes  the 
TTM  equation  for  the  lattice  temperature.  The  diffusion 
equation  for  the  electron  temperature  is  solved  by  a  finite- 
difference  method  simultaneously  with  MD  integration  of 
the  equations  of  motion  of  atoms.  The  electron  temperature 
enters  a  coupling  term  that  is  added  to  the  MD  equations  of 
motion  to  account  for  the  energy  exchange  between  the  elec¬ 
trons  and  the  lattice.  The  cells  in  the  finite-difference  dis¬ 
cretization  are  related  to  the  corresponding  volumes  of  the 
MD  system  and  the  local  lattice  temperature  is  defined  for 
each  cell  from  the  average  kinetic  energy  of  thermal  motion 
of  atoms.  A  complete  description  of  the  combined  TTM-MD 
model  is  given  in  Ref.  19. 


Irradiation  by  a  laser  pulse  is  represented  in  the  con¬ 
tinuum  part  of  the  model  by  a  source  term  with  a  Gaussian 
temporal  profile  and  exponential  attenuation  of  laser  inten¬ 
sity  with  depth  under  the  surface  (Beer-Lambert  law).  The 
electron  mean  free  path  in  Au  is  larger  than  the  optical  pen¬ 
etration  depth  and  the  ballistic  energy  transport  defines  the 
effective  laser  energy  deposition  depth,  estimated  to  be  on 
the  order  of  100  nm.26  Since  the  ballistic  range  in  Au  ex¬ 
ceeds  the  thickness  of  the  films  considered  in  this  work, 
20  nm,  the  reflection  of  the  ballistic  electrons  from  the  back 
surface  of  the  film  results  in  a  uniform  distribution  of  the 
electronic  temperature  established  on  the  time  scale  of  elec¬ 
tron  thermalization.  The  theoretical  prediction  of  the  uniform 
rise  of  the  electronic  temperature  for  thicknesses  smaller  than 
the  ballistic  range  has  been  confirmed  in  a  series  of  pump- 
probe  measurements  of  transient  reflectivity  performed  for 
Au  films  of  different  thicknesses,  from  10  to  500  nm.27  The 
effect  of  the  ballistic  energy  transport  and  the  finite  size  of 
the  film  are  accounted  for  in  the  source  term  describing  the 
laser  irradiation.19 

The  range  of  laser  fluences  used  in  the  simulations  per¬ 
formed  for  Au  films,  from  45  to  180  J/m2,  is  chosen  so  that 
only  an  incomplete  heterogeneous  melting  of  the  film  is  ob¬ 
served  at  the  lowest  fluence  and  an  ultrafast  homogeneous 
melting  of  the  whole  film  is  observed  at  the  highest  fluence. 
The  absorbed  laser  fluences  rather  than  the  incident  fluences 
are  given  here  and  are  used  in  the  remaining  part  of  the 
paper. 

The  interatomic  interaction  in  the  MD  part  of  the  model  is 
described  by  the  embedded-atom  method  (EAM)  with  the 
functional  form  and  parametrization  suggested  in  Ref.  28. 
The  choice  of  the  interatomic  potential  defines  all  the  ther¬ 
mal  and  elastic  properties  of  the  material.  Some  of  the  prop¬ 
erties  of  EAM  Au  relevant  to  the  material  response  to  the 
laser  heating  are  listed  in  Table  I,  along  with  experimental 
data  for  Au.  While  there  are  some  quantitative  discrepancies 
between  the  properties  of  the  model  EAM  Au  and  experi¬ 
mental  data,  the  overall  agreement  is  reasonable  and  we  can 
expect  that  the  model  will  adequately  reproduce  the  material 
response  to  fast  laser  heating.  Moreover,  the  knowledge  of 
the  thermodynamic  parameters  of  the  model  material  allows 
us  to  perform  a  quantitative  analysis  and  physical  interpreta¬ 
tion  of  the  simulation  results. 

The  parameters  used  in  the  TTM  equation  for  the  electron 
temperature  are  as  follows.  For  Au,19  the  electronic  heat  ca¬ 
pacity  is  Ce=yTe  with  y=  7 1  J  m-3  K-2,  the  electron-phonon 
coupling  constant  is  G  =  2.1  X  1016  W  m-3  K-1,  and  the  de¬ 
pendence  of  the  electron  thermal  conductivity  on  the  electron 
and  lattice  temperatures  is  described  by  an  expression  sug¬ 
gested  in  Ref.  32.  For  Al,  Ce=yTe  with  y=  1 25  J  m~3  K-2, 
G  =  3.1  X  1017  W  m-3  K-1,33  the  electron  thermal  conductiv¬ 
ity  is  Ke  =  KQTe/Ti  with  A0  =  238  W  m-  K-1,34  and  the  opti¬ 
cal  penetration  depth  at  800  nm  is  8  nm.12 

The  initial  MD  system  in  simulations  of  laser  interaction 
with  Au  films  is  a  fee  crystal  composed  of  500  000  atoms 
with  dimensions  of  20.46  X  20.46  X  20.46  nm3  and  periodic 
boundary  conditions  imposed  in  the  directions  parallel  to  two 
(001)  free  surfaces.  The  periodic  boundary  conditions  simu¬ 
late  the  situation  in  which  the  laser  spot  diameter  is  suffi¬ 
ciently  large  so  that  the  energy  redistribution  in  the  lateral 
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TABLE  I.  Some  of  the  material  parameters  determined  for  the  EAM  Au  material.  Values  of  the  equilib¬ 
rium  melting  temperature  Tm,  volume  change  AT,,,,  enthalpy  A Hm,  and  entropy  A5„,  of  melting  are  given  for 
zero  pressure.  The  dependence  of  the  equilibrium  melting  temperature  on  pressure  ( dT/dP)m  is  determined 
from  liquid-crystal  coexistence  simulations  and  confirmed  by  the  calculations  based  on  the  Clapeyron  equa¬ 
tion  (dT/dP)m=Wml ASm.  The  value  given  in  the  table  is  calculated  for  zero  pressure.  Variations  of  the 
coefficient  of  linear  expansion,  a,  and  heat  capacity  at  zero  pressure,  Cp,  are  given  for  a  temperature  range 
from  293  to  950  K.  Experimental  values  for  Au  are  from  Refs.  29-31. 


An 

(K) 

AV,„ 

(cm3  mol  *) 

A  Sm 

(J  K_1  mol-1) 

A  Hm 

(kj  mol  *) 

(dT/dP)m 
(K  GPa-1) 

Cp 

(J  K_1  mol-1) 

a 

(10“6  K”1) 

EAM  Au 

963 

0.28 

8.7 

8.4 

32.2 

25.9-30.3 

10.3-21.4 

Experiment 

1336 

0.55 

9.6 

12.8 

57.5 

25.4-31.2 

14.2-19.1 

Ref.  29 

Ref.  29 

Ref.  29 

Ref.  29 

Ref.  29 

Ref.  30 

Ref.  31 

directions,  parallel  to  the  free  surfaces  of  the  film,  can  be 
neglected  on  the  time  scales  considered  in  the  simulations. 
Several  simulations  are  performed  for  systems  with  dimen¬ 
sions  of  8.18  X  8.18  X  20.46  nm3  (80  000  atoms),  16.37 
X  6.37  X  20.46  nm3  (320  000  atoms),  and  28.64X28.64 
X  20.46  nm3  (980  000  atoms)  to  investigate  the  effect  of  the 
size  of  the  computational  cell  on  the  calculated  diffraction 
profiles.  In  the  simulation  performed  for  an  A1  film,  a  similar 
system  composed  of  500  000  atoms  with  dimensions  of 
20.56  X  20.56  X  20.56  nm3  is  used.  Before  applying  laser  ir¬ 
radiation,  all  systems  are  equilibrated  at  300  K  and  zero 
pressure. 

The  identification  of  liquid  and  crystal  regions  in  the 
atomic  configurations  obtained  in  the  simulations  is  done 
with  a  local  order  parameter  calculated  for  each  atom  based 
on  the  local  structure  within  the  first  two  neighbor  shells.19 
The  local  order  parameter  is  used  to  identify  the  crystal  and 
liquid  regions  in  the  transient  atomic  configurations  and  to 
quantitatively  describe  the  kinetics  of  the  melting  process. 


III.  NUMERICAL  METHODS  FOR  CALCULATION  OF 
STRUCTURE  FUNCTIONS  FROM  ATOMIC 
CONFIGURATIONS 

In  order  to  calculate  diffraction  patterns  from  atomic  con¬ 
figurations  generated  in  MD  simulations,  we  consider  the 
scattering  of  a  monochromatic  or  monoenergetic  beam  of 
x-ray  photons  or  electrons  on  a  sample  consisting  of  N  at¬ 
oms.  Assuming  that  only  single  elastic  scattering  takes  place, 
the  amplitude  of  the  wave  scattered  by  the  sample  is  given 
by  summing  the  amplitudes  of  scattering  from  each  atom  in 
the  configuration:35,36 

N 

( 2)  =  2  fi  exp (-  <  2  •  r,) ,  ( 1 ) 

i=  l 

where  Q  is  the  scattering  vector,  r,  is  the  position  of  atom  i 
with  respect  to  an  arbitrarily  chosen  origin,  /,  is  the  x-ray  or 
electron  atomic  scattering  form  factor  for  the  z'th  atom,  and 
the  sum  goes  over  all  the  atoms.  The  magnitude  of  Q  is  given 
by  <2=477  sin  6/  X,  where  6  is  half  the  angle  between  the 
incident  and  scattered  wave  vectors,  and  X  is  the  wavelength 
of  the  incident  wave.  The  scattering  form  factors  /,  are  func¬ 


tions  of  Q  with  different  dependences  in  x-ray  and  electron 
scattering.37  The  intensity  of  the  scattered  wave  can  be  found 
by  multiplying  the  scattered  wave  function  by  its  complex 
conjugate: 

N  N 

AQ)  =  ^(2)  •  ^*(2)  =  2  2  fifj  exp[-  iQ  ■  (n -  fj)]. 

7=1  i=  i 

(2) 


The  spherically  averaged  powder-diffraction  intensity  profile 
can  be  obtained  by  integration  of  Eq.  (2)  over  all  directions 
of  interatomic  separation  vector  rij=ri-rj,  resulting  in  the 
Debye  scattering  equation35 


N  N 


/(2)  =  22/i/> 

7=1  <=  1 


sin(2u;) 

Qrij 


(3) 


By  dividing  the  above  equation  by  we  obtain  a  func¬ 

tion  that,  following  Ref.  38,  we  call  the  structure  function: 


5(2)  = 


AQ)  _ .  -  y  y  ff  sin(2u;) 

2/-  2/fli<7' 

i=l  /=1 


(4) 


This  function  approaches  unity  at  large  Q  and  is  often  used 
to  present  experimental  diffraction  results.  For  a  monatomic 
system  the  dependence  on  the  atomic  form  factors  can  be 
eliminated  and  we  have 


5(2)  =  1  +  fr2  2 


N 


7=1  i<7 


sin(2u;) 

Qrij 


(5) 


The  structure  function  defined  by  Eq.  (5)  can  be  computed 
directly  from  the  atomic  configurations  generated  in  MD 
simulations.  The  calculations,  however,  involve  the  summa¬ 
tion  over  all  pairs  of  atoms  in  the  system,  leading  to  the 
quadratic  dependence  of  the  computational  cost  on  the  num¬ 
ber  of  atoms  and  making  the  calculations  prohibitively  ex¬ 
pensive  for  large  systems.39 

An  alternative  approach  to  calculation  of  5(2)  is  to  sub¬ 
stitute  the  double  summation  over  atomic  positions  in  Eq.  (5) 
by  integration  over  the  pair  density  function,  which  is  a  real- 
space  representation  of  correlations  in  atomic  positions,38 
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where  8  is  the  Dirac  delta  function.  Although  the  calculation 
of  the  pair  density  function  still  involves  N2/  2  evaluations  of 
interatomic  distances  rVp  it  can  be  done  much  more  effi¬ 
ciently  than  the  double  summation  in  Eq.  (5),  which  requires 
evaluation  of  the  sine  function  and  repetitive  calculations  for 
each  value  of  Q.  The  expression  for  the  structure  function, 
Eq.  (5),  can  be  now  reduced  to  a  simple  integration,  which, 
in  fact,  is  the  Fourier  transform  of  the  pair  density  function: 

f"  -i  sin(Qr) 

S(Q)  =  1+  4t TrpW—^dr.  (7) 

Jo  Qr 


In  the  calculation  of  the  pair  density  function,  the  maxi¬ 
mum  value  of  r  is  limited  by  the  size  of  the  MD  computa¬ 
tional  cell.  If  the  periodic  boundary  conditions  are  used  to 
represent  a  part  of  a  larger  or  infinite  system,  the  pair  density 
function  can  only  be  evaluated  up  to  a  maximum  value  /?max 
that  should  not  exceed  one-half  of  the  computational  cell. 
The  truncation  of  the  numerical  integration  in  Eq.  (7)  at  7?max 
induces  spurious  ripples  with  a  period  of  A  =  2tt/ Rmax.40  A 
number  of  methods  have  been  proposed  to  suppress  these 
ripples  so  that  the  Fourier  ringing  dies  out  more 
quickly. 35,38’39,41'42  The  method  that  we  adopt  in  this  work  is 
to  multiply  the  integrand  in  Eq.  (7)  by  a  damping  function 
W{r),  similar  to  the  Lorch  modification  function  in  neutron 
diffraction  experiments,43 


W(r)  = 


rA 


2 


3in(v~) 


Rn 


(8) 


Thus  the  structure  function  S(Q )  can  be  calculated  as 


S{Q )  =  1  + 


l 


„  7  ^  sln(Gr) 

4ttt  p(r) — — - 

o  Qr 


W(r)dr. 


(9) 


The  damping  function  replaces  the  sharp  step  function  at  the 
cutoff  distance  7?max  by  a  smoothly  decreasing  contribution 
from  the  density  function  at  large  interatomic  distances  and 
eventually  approaching  zero  at  7?max. 

Figure  1(a)  shows  structure  functions  calculated  with  and 
without  the  damping  function  for  an  fee  Au  system  com¬ 
posed  of  500  000  atoms  (/?max=  100  A)  and  equilibrated  at 
300  K.  It  is  apparent  that  the  introduction  of  the  damping 
function  completely  eliminates  spurious  truncation  ripples  in 
the  structure  function.  While  for  the  fee  crystallite  at  300  K 
the  presence  of  ripples  does  not  prevent  identification  of  the 
real  structural  peaks,  Fig.  1(a),  the  elimination  of  ripples  is 
crucial  for  the  analysis  of  the  structural  transformations  oc¬ 
curring  at  elevated  temperatures,  when  the  real  structural 
peaks  can  be  small  and  completely  obscured  by  the  ripples. 

Although  application  of  the  damping  function  eliminates 
the  truncation  ripples,  the  finite  size  of  the  MD  system  and 
the  introduction  of  the  cutoff  in  the  pair  density  function  also 


(b)  Q  (A  ’) 


FIG.  1.  (Color  online)  The  effect  of  the  finite  size  of  the  MD 
system  on  the  calculated  structure  functions.  Elimination  of  spuri¬ 
ous  ripples  induced  by  the  truncation  of  the  pair  density  function  at 
Rmax=100  A  by  introduction  of  the  damping  function  W(r)  in  Eq. 
(9)  is  illustrated  in  (a),  where  the  results  are  shown  for  a  20.46 
X  20.46  X  20.46  nm3  fee  Au  system  equilibrated  at  300  K.  Solid 
(black)  and  dashed  (red)  lines  show  the  results  of  the  calculations 
performed  with  and  without  W(r),  respectively.  Structure  functions 
calculated  using  Eq.  (9),  with  the  damping  function,  for  systems  of 
four  different  sizes  in  the  lateral  directions  are  shown  in  (b). 

affect  the  real  peaks  of  the  structure  function.  To  illustrate 
the  effect  of  the  size  of  the  system  on  the  characteristics  of 
structure  function,  the  results  of  calculations  performed  for  a 
20-nm-thick  Au  film  represented  by  MD  computational  cells 
with  four  different  sizes  in  the  directions  of  periodic  bound¬ 
ary  conditions,  parallel  to  the  surfaces  of  the  film,  8.18, 
16.37,  20.46,  and  28.64  nm  are  shown  in  Fig.  1(b).  The  cut¬ 
off  distances  Rmax,  used  in  the  calculations  of  the  pair  density 
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functions  are  40,  80,  100,  and  140  A,  respectively.  It  can  be 
seen  that  the  failure  to  include  long-range  atomic  correla¬ 
tions  beyond  40  A  results  in  a  significant  broadening  of  all 
peaks,  with  some  of  the  peaks  starting  to  merge.  As  the  value 
of  Rmax  increases  to  80  A,  the  peaks  become  sharp  and  well 
defined.  Further  increase  of  /?max  to  100  and  to  140  A  results 
in  a  much  more  moderate  sharpening  of  the  peaks,  with  no 
changes  to  the  peak  positions.  Analytical  calculation  of  the 
broadening  of  the  peaks  due  to  the  introduction  of  the  damp¬ 
ing  function  with  R„mx  =  1 00  A  predicts  a  broadening  of 
A(9  =  5.437//?max=0.054  37  A-1.44  Since  the  purpose  of  the 
present  paper  is  to  investigate  the  evolution  of  the  diffraction 
pattern  during  the  fast  laser  heating  of  the  film  up  to  the 
melting  temperature  and  above,  the  relatively  small  finite- 
size  effect  is  expected  to  be  negligible  compared  to  the 
changes  in  the  diffraction  peaks  associated  with  the  tempera¬ 
ture  rise  and  structural  transformations.  Thus,  all  simulations 
discussed  in  the  next  section  are  performed  for  a  20.46 
X  20.46  X  20.46  nm3  (500  000  atoms)  system  and  all  struc¬ 
ture  functions  are  calculated  with  Eq.  (9)  and  Ajnax  =  100  A. 


IV.  RESULTS  AND  DISCUSSION 
A.  Kinetics  and  mechanisms  of  laser  melting 

The  time  scales  of  laser  melting  predicted  in  TTM-MD 
simulations  of  20  nm  Au  films  irradiated  with  200  fs  laser 
pulses  at  absorbed  fluences  ranging  from  45  to  180  J/m2  are 
presented  in  Fig.  2.  The  fraction  of  the  crystal  phase  is  de¬ 
fined  by  the  number  of  atoms  with  local  crystalline  environ¬ 
ment,  as  predicted  by  the  local  order  parameter.19  Two  dis¬ 
tinct  regimes  can  be  identified  in  Fig.  2,  a  high-fluence 
regime  when  the  entire  film  melts  within  just  several  pico¬ 
seconds,  and  a  low-fluence  regime  when  the  melting  process 
slows  down  and  the  melting  starting  time  increases  sharply 
with  decreasing  fluence.  Below  we  briefly  discuss  the  melt¬ 
ing  mechanisms  of  Au  films  in  these  two  regimes. 

At  the  highest  laser  fluence  of  180  J/m2,  melting  starts  at 
about  10  ps  after  the  laser  pulse  and  completes  by  14  ps. 
Similarly  fast  decrease  of  the  fraction  of  the  crystal  phase  is 
observed  in  simulations  performed  at  absorbed  fluences  of 
140  and  100  J/m2.  The  melting  mechanism  in  this  high- 
fluence  regime  is  exemplified  here  by  the  results  obtained  in 
a  simulation  performed  at  180  J/m2.  The  temporal  and  spa¬ 
tial  evolution  of  the  lattice  temperature  and  pressure  in  the 
irradiated  film  is  shown  in  the  form  of  contour  plots  in  Fig. 
3.  The  large  mean  free  path  of  excited  electrons  and  a  weak 
electron-phonon  coupling  in  Au  results  in  a  uniform  distri¬ 
bution  of  the  electronic  temperature  in  the  film  before  the 
electron-lattice  thermalization.  As  a  result  the  whole  film  is 
heated  up  uniformly  at  a  rate  on  the  order  of  ~  1014  K/ s  [Fig. 
3(a)],  The  fastest  heating  is  observed  within  the  first  15  ps  of 
the  simulation,  whereas  complete  equilibration  between  the 
hot  electrons  and  the  lattice  takes  up  to  70  ps.  The  fast  in¬ 
crease  of  the  lattice  temperature  results  in  the  overheating  of 
the  lattice  above  the  limit  of  its  stability,  ~1.257j„,19  and 
leads  to  the  fast  homogeneous  melting  of  the  whole  film 
within  ~3-4  ps.  The  solid  and  dashed  lines  in  Fig.  3  mark 
the  beginning  (90%  of  the  crystal  phase)  and  the  end  (10%  of 


FIG.  2.  (Color  online)  The  time  scales  of  the  melting  process  in 
a  20  nm  Au  film  irradiated  with  a  200  fs  laser  pulse  at  different 
absorbed  fluences.  Semilogarithmic  plots  of  the  time  required  to 
melt  certain  fractions  of  the  film  are  shown  in  (a).  Each  curve 
corresponds  to  a  particular  fraction  of  the  remaining  crystal  phase 
as  a  function  of  the  absorbed  fluence.  For  example,  the  (red)  curve 
with  squares  corresponds  to  the  time  after  laser  excitation  when 
90%  of  the  atoms  in  the  film  belong  to  the  crystal  phase.  The  de¬ 
crease  of  the  fraction  of  the  crystal  phase  with  time  is  shown  for 
each  simulation  in  (b).  The  atoms  in  the  crystal  phase  are  distin¬ 
guished  from  the  ones  in  the  liquid  phase  based  on  the  local  order 
parameter  (Ref.  19). 

the  crystal  phase)  of  the  melting  process  and  can  be  related 
to  the  corresponding  points  in  Fig.  2. 

The  pressure  contour  plot  in  Fig.  3(b)  shows  that  the  fast 
lattice  heating  results  in  the  buildup  of  compressive  stresses 
inside  the  film  within  the  first  ~5  ps.  The  initial  compressive 
pressure  drives  the  expansion  of  the  film.  In  earlier  simula¬ 
tions  performed  for  Ni  and  thicker  Au  films19,20,45  the  relax  - 
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FIG.  3.  (Color  online)  Contour  plots  of  the  lattice  temperature 
(a)  and  pressure  (b)  for  a  simulation  of  laser  melting  of  a  20  nrn  Au 
film  irradiated  with  a  200  fs  laser  pulse  at  an  absorbed  fluence  of 
180  J/m2.  Solid  and  dashed  lines  show  the  beginning  and  the  end 
of  the  melting  process.  Laser  pulse  is  directed  along  the  Y  axes, 
from  the  top  of  the  contour  plots.  The  stepwise  shape  of  the  contour 
plot  boundaries  is  related  to  the  discretization  of  the  mesh  over 
which  average  temperature  and  pressure  values  are  calculated. 

ation  of  the  laser-induced  compressive  stresses  resulted  in 
generation  of  tensile  stresses,  followed  by  pressure  oscilla¬ 
tions  or  even  disintegration  of  the  film.  For  a  20  nm  free¬ 
standing  Au  film,  however,  the  time  of  the  mechanical  relax¬ 
ation  is  on  the  order  of  5  ps  (time  needed  for  two  unloading 
waves  propagating  from  the  free  surfaces  of  the  film  to  cross 
one  half  of  the  depth  of  the  film),  significantly  shorter  than 
the  time  of  the  lattice  heating.  As  a  result,  the  film  expands 
during  the  lattice  heating  and  pressure  (and  film  thickness) 
oscillations  can  be  hardly  observed  in  Fig.  3(b). 

The  atomic-level  picture  of  the  homogeneous  melting 
process  in  the  high-energy  regime  is  shown  in  Fig.  4(a).  The 
visual  analysis  of  snapshots  from  the  simulation  suggests 
that  the  growth  of  the  liquid  regions  appearing  at  the  free 


surfaces  of  the  film  (see  a  snapshot  taken  at  10  ps)  does  not 
make  any  significant  contribution  to  the  overall  melting  pro¬ 
cess.  The  energy  transfer  from  the  hot  electrons  to  the  lattice 
quickly  leads  to  the  overheating  of  the  lattice  up  to  the  limit 
of  the  crystal  stability,  when  a  spontaneous  nucleation  of  a 
large  number  of  small  liquid  regions  occurs  throughout  the 
film,  leading  to  the  rapid  collapse  of  the  crystalline  structure 
from  10  to  14  ps. 

A  similar  melting  process,  dominated  by  homogeneous 
nucleation  of  liquid  regions  inside  the  overheated  crystal,  is 
observed  for  all  other  fluences  in  the  high-fluence  regime, 
down  to  70  J/m2.  The  decrease  in  fluence  shifts  the  melting 
process  to  a  later  time  and  gradually  increases  the  contribu¬ 
tion  of  the  heterogeneous  melting  that  takes  place  by  propa¬ 
gation  of  two  melting  fronts  from  the  surfaces  of  the  film. 
The  latter  process  is  reflected  in  the  development  of  the  ini¬ 
tial  slow  shoulders  in  the  melting  curves  shown  in  Fig.  2(b). 
The  slow  components  of  the  melting  process  are  apparent  in 
the  curve  plotted  for  70  J/m2,  where  almost  10%  of  the  film 
melts  due  to  the  propagation  of  the  melting  fronts,  before  the 
critical  superheating  is  reached  and  the  faster  homogeneous 
melting  takes  over. 

The  melting  process  observed  at  55  J/m2  can  be  consid¬ 
ered  to  be  a  transitional  one  between  the  high-fluence  (ho¬ 
mogeneous  melting)  and  low-fluence  (heterogeneous  melt¬ 
ing)  regimes.  At  this  fluence  both  mechanisms  of  melting 
contribute  approximately  equally  to  the  melting  process  [Fig. 
5(a)].  During  the  time  from  25  to  40  ps,  the  propagation  of 
two  melting  fronts  from  the  surfaces  of  the  film  is  largely 
responsible  for  the  increase  in  the  fraction  of  the  liquid 
phase.  After  40  ps,  the  growth  of  liquid  regions  nucleated 
inside  the  overheated  crystal  starts  to  make  a  major  contri¬ 
bution  to  the  melting  process,  significantly  accelerating  the 
total  rate  of  melting,  as  reflected  in  the  slope  of  the  melting 
line  in  Fig.  2(b).  The  melting  slows  down  again  at  later  time 
due  to  the  lattice  temperature  decrease.  From  the  time  of  70 
to  90  ps  the  last  crystalline  island  slowly  melts  [Fig.  5(a)],  as 
the  temperature  of  the  film  is  approaching  the  equilibrium 
melting  temperature. 

At  laser  fluences  below  55  J/m2,  the  propagation  of  the 
melting  fronts  from  the  surfaces  of  the  film  dominates  the 
melting  process  and  a  significant  increase  of  the  melting  time 
is  observed  in  Fig.  2.  Snapshots  from  a  simulation  performed 
at  a  fluence  of  45  J/m2  illustrate  the  melting  process  in  this 
regime  [Fig.  6(a)].  Two  melting  fronts  propagate  from  the 
free  surfaces  with  velocities  that  decrease  down  to  zero  as 
the  melting  progresses  and  the  temperature  approaches  the 
equilibrium  melting  temperature.  By  the  end  of  the  simula¬ 
tion  (500  ps)  the  temperature  of  the  film  decreases  down  to 
the  equilibrium  melting  temperature  and  the  remaining  49% 
of  the  film  material  remains  in  the  crystalline  state.  Incom¬ 
plete  and  purely  heterogeneous  melting  is  also  observed  in 
the  simulations  performed  at  50  and  51  J/m2.  Based  on  the 
properties  of  the  model  EAM  Au  material  we  can  estimate 
the  critical  absorbed  fluence  Fm  that  would  supply  enough 
energy  to  completely  melt  the  film: 

Fm  =(|  C,(T)dT+  f  Ce(T)dT+AHm)d,  (10) 

\  J  300  J  300  / 

where  d  is  the  thickness  of  the  film,  Q  and  Ce  are  the  lattice 
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FIG.  4.  (Color  online)  Snapshots  (a)  and 
structure  functions  (b)  of  atomic  configurations 
during  the  melting  process  in  a  20  nm  Au  film 
irradiated  with  a  200  fs  laser  pulse  at  an  absorbed 
fluence  of  180  J/m2.  Atoms  are  colored  accord¬ 
ing  to  the  local  order  parameter  </> — blue  atoms 
have  local  crystalline  surroundings,  red  atoms  be¬ 
long  to  the  liquid  phase.  In  (a),  the  laser  pulse  is 
directed  from  the  right  to  the  left  sides  of  the 
snapshots.  In  (b),  curves  are  shifted  vertically 
with  respect  to  each  other  in  order  to  better  show 
the  changes  in  the  structure  function.  Zero  time 
corresponds  to  a  perfect  fee  crystal  at  300  K  just 
before  the  laser  irradiation. 


and  electron  heat  capacities,  and  A Hm  is  the  heat  of  melting. 
Using  the  parameters  of  EAM  Au  listed  in  Table  I  and  elec¬ 
tron  heat  capacity  defined  in  Sec.  II  we  can  estimate  Fm 
=  53.78  J/m2,  with  36.52  and  0.59  J/m2  going  into  heating 
of  the  lattice  and  electrons  up  to  the  equilibrium  melting 


temperature  Tm  and  16.67  J/m2  required  for  melting  the  film 
at  Tm.  This  estimation  is  consistent  with  the  results  of  MD 
simulation  described  above. 

The  relative  contribution  of  the  homogeneous  and  hetero¬ 
geneous  melting  mechanisms  in  the  simulations  described 
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FIG.  5.  (Color  online)  Snapshots  (a)  and 
structure  functions  (b)  of  atomic  configurations 
during  the  melting  process  in  a  20  nm  Au  film 
irradiated  with  a  200  fs  laser  pulse  at  an  absorbed 
fluence  of  55  J/m2.  Atoms  are  colored  according 
to  the  local  order  parameter  <f> — blue  atoms  have 
local  crystalline  surroundings,  red  atoms  belong 
to  the  liquid  phase.  In  (a),  the  laser  pulse  is  di¬ 
rected  from  the  right  to  the  left  sides  of  the  snap¬ 
shots.  In  (b),  curves  are  shifted  vertically  with 
respect  to  each  other  in  order  to  better  show  the 
changes  in  the  structure  function.  Zero  time  cor¬ 
responds  to  a  perfect  fee  crystal  at  300  K  just 
before  the  laser  irradiation. 


above  is  controlled  by  the  temperature  dependence  of  the 
velocity  of  the  melting  fronts  propagating  from  the  free  sur¬ 
faces  of  the  film  and  the  lattice  heating  regime.  The  rate  of 
the  lattice  heating  is  defined  by  the  irradiation  parameters 
(laser  fluence  and  pulse  duration)  and  the  strength  of  the 


electron-phonon  coupling,  whereas  the  temperature  depen¬ 
dence  of  the  velocity  of  the  melting  front  can  be  described  by 
nonequilibrium  kinetic  theory.46  The  maximum  velocity  of 
the  melting  front  propagation  is  the  velocity  at  the  limit  of 
the  crystal  stability,  above  which  a  massive  homogeneous 
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FIG.  6.  (Color  online)  Snap¬ 
shots  (a)  and  structure  functions 
(b)  of  atomic  configurations  dur¬ 
ing  the  melting  process  in  a  20  nm 
Au  film  irradiated  with  a  200  fs 
laser  pulse  at  an  absorbed  fluence 
of  45  J/m2.  Atoms  are  colored  ac¬ 
cording  to  the  local  order  param¬ 
eter  <f) — blue  atoms  have  local 
crystalline  surroundings,  red  at¬ 
oms  belong  to  the  liquid  phase.  In 
(a),  the  laser  pulse  is  directed 
from  the  right  to  the  left  sides  of 
the  snapshots.  In  (b),  curves  are 
shifted  vertically  with  respect  to 
each  other  in  order  to  better  show 
the  changes  in  the  structure 
function. 


nucleation  of  liquid  regions  inside  the  overheated  crystal 
takes  place.  While  one  can  safely  assume  that  the  maximum 
velocity  of  the  melting  front  is  ultimately  limited  by  the 
speed  of  sound,47  recent  MD  simulations19  demonstrate  that 
the  real  maximum  velocity  at  the  limit  of  crystal  stability 
does  not  exceed  15%  of  the  speed  of  sound.  Moreover,  under 
conditions  of  laser  melting,  the  overheating  required  for  the 


onset  of  the  homogeneous  nucleation  of  liquid  regions  is 
significantly  reduced  by  the  uniaxial  lattice  distortions  pro¬ 
duced  as  a  result  of  the  relaxation  of  laser-induced  ther¬ 
moelastic  stresses.20  This  reduction  in  the  lattice  stability  ex¬ 
plains  why  the  homogeneous  nucleation  is  observed  down  to 
a  relatively  low  fluence  of  55  J/m2,  when  the  melting  pro¬ 
cess  is  slow  and  takes  up  to  90  ps. 
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B.  Structure  function  calculations 

Calculation  of  the  diffraction  profiles  provides  an  attrac¬ 
tive  possibility  to  directly  relate  the  detailed  information  on 
the  kinetics  and  mechanisms  of  laser-induced  structural 
transformations  obtained  in  MD  simulations  to  the  results  of 
time-resolved  diffraction  experiments. 5-8,10-13  In  this  section 
we  present  the  results  on  the  evolution  of  the  structure  func¬ 
tion  S{Q)  in  two  representative  simulations  discussed  above, 
a  simulation  at  1 80  J  /  m2  where  a  fast  homogeneous  melting 
is  observed  and  a  simulation  at  45  J  /  nr  where  much  slower 
heterogeneous  melting  takes  place. 

The  temporal  evolution  of  the  diffraction  profile  is  shown 
for  a  fluence  of  180  J/m2  in  Fig.  4(b).  Although  by  the  time 
of  7  ps  the  melting  process  has  barely  started  and  most  at¬ 
oms  (—99%)  are  still  identified  by  the  local  order  parameter 
as  maintaining  local  crystalline  surroundings,  the  structure 
function  has  changed  significantly.  First,  there  is  a  consider¬ 
able  reduction  in  the  heights  of  all  the  peaks  as  compared  to 
the  structure  function  calculated  at  zero  time,  before  the  laser 
pulse.  This  reduction  can  be  largely  attributed  to  the  fast 
heating  of  the  crystal  from  300  to  1059  K  [Fig.  3(a)].  Quan¬ 
titative  analysis  of  peak  height  reduction  due  to  the  increase 
in  the  amplitudes  of  thermal  atomic  vibrations  is  given  in 
Sec.  IV  D.  Second,  there  is  a  pronounced  shift  to  the  left  of 
the  (111)  diffraction  peak  and  splitting  of  the  (200),  (220), 
and  (311)  peaks.  The  appearance  of  new  diffraction  peaks 
may  be  indicative  of  solid-solid  phase  transformations.  A 
detailed  structural  analysis,  however,  does  not  reveal  any 
structural  changes  in  the  film  beyond  the  appearance  of  a 
relatively  small  number  of  point  defects  (vacancy-interstitial 
pairs).  Both  the  shift  and  the  splitting  of  the  peaks  are  actu¬ 
ally  related  to  the  uniaxial  expansion  of  the  film  in  response 
to  the  laser  heating,  as  explained  in  Sec.  IV  C. 

A  fast  collapse  of  the  crystalline  structure  from  10  to 
14  ps,  apparent  from  the  snapshots  shown  in  Fig.  4(a),  is  also 
clearly  reflected  in  the  evolution  of  the  structure  function.  In 
just  3-4  ps  the  peaks  characteristic  of  the  fee  structure  dis¬ 
appear  and  the  structure  function  takes  the  shape  character¬ 
istic  of  the  liquid  structure,  with  only  one  broad  peak  that 
can  be  identified.  The  duration  of  the  melting  process  ob¬ 
served  in  this  simulation  is  in  a  good  agreement  with  the  one 
measured  in  time-resolved  electron  diffraction  experiments 
performed  for  20  nm  A1  films,  3.5  ps.11  Some  of  the  differ¬ 
ences  between  the  evolution  of  the  diffraction  peaks  ob¬ 
served  in  the  simulation  and  experiment,  most  notably  the 
absence  of  peak  splitting  and  shifts  in  the  positions  of  dif¬ 
fraction  peaks  in  experimental  observations,  are  discussed  in 
Sec.  V. 

An  evolution  of  the  structure  function  in  the  low-fluence 
regime,  when  the  melting  process  is  dominated  by  the  propa¬ 
gation  of  the  melting  fronts  from  the  free  surfaces  of  the 
film,  is  illustrated  in  Fig.  6(b),  where  the  results  are  shown 
for  an  absorbed  fluence  of  45  J/m2.  Similarly  to  the  fast 
homogeneous  melting  discussed  above,  laser  irradiation 
leads  to  the  reduction  of  heights  of  the  peaks  and  induces  a 
shift  of  the  (111)  peak  to  the  left  and  splitting  of  other  peaks. 
There  are,  however,  important  differences  in  the  evolution  of 
the  structure  functions  in  low-  and  high-energy  regimes. 
While  the  heights  of  the  peaks  are  decreasing  gradually  dur- 


TABLE  II.  Fraction  of  atoms  with  local  crystalline  surroundings 
(local  order  parameter  </>>0.04)  in  a  20  nm  Au  film  irradiated  with 
a  200  fs  laser  pulse  at  absorbed  fluences  of  45  and  180  J/m2. 


Absorbed  fluence  =  45  J/m2 


Time  (ps) 

10 

60 

180 

500 

Crystal  fraction  (%) 

100.0 

93.1 

67.5 

49.1 

Absorbed  fluence = 

180  J/m2 

Time  (ps) 

11 

12 

13 

14 

Crystal  fraction  (%) 

93.1 

67.5 

15.9 

0.3 

ing  the  melting  process,  the  peaks  remain  well  defined  up  to 
the  end  of  the  simulation,  when  more  than  50%  of  the  film  is 
melted.  For  the  same  fraction  of  atoms  belonging  to  the  crys¬ 
talline  parts  of  the  system  (having  crystalline  local  environ¬ 
ment  as  defined  by  the  local  order  parameter)  the  diffraction 
peaks  are  more  pronounced  in  the  case  of  heterogeneous 
melting.  For  instance,  although  from  Table  II  we  can  see  that 
the  same  number  of  atoms  remain  in  the  crystalline  parts  of 
the  film  at  11  ps  after  the  180  J/m2  pulse  and  at  60  ps  after 
the  45  J/m2  pulse,  the  diffraction  peaks  are  much  more 
clearly  defined  in  the  corresponding  curve  in  Fig.  6(b),  as 
compared  to  the  one  in  Fig.  4(b).  One  can  derive  the  same 
conclusion  from  the  comparison  of  structure  functions  shown 
for  12  ps  in  Fig.  4(b)  and  180  ps  in  Fig.  6(b). 

It  is  clear  from  the  visual  analysis  of  the  snapshots  of 
atomic  configurations  that  the  homogeneous  nucleation  of  a 
large  number  of  liquid  regions  throughout  the  film  [Fig. 
4(a)],  is  effective  in  destabilizing  the  lattice  and  reducing  the 
long-range  order  throughout  the  film.  While  most  of  the  at¬ 
oms  still  retain  their  local  crystalline  surroundings  at  1 1  and 
12  ps,  the  long-range  correlations  in  atomic  positions  are 
largely  lost  and  the  diffraction  peaks  are  significantly  re¬ 
duced.  Similarly,  the  presence  of  relatively  small  crystalline 
islands  observed  at  60  and  70  ps  in  a  simulation  performed  at 
55  J/m2  [Fig.  5(a)]  can  be  hardly  identified  from  the  corre¬ 
sponding  diffraction  patterns  [Fig.  5(b)].  In  the  case  of  het¬ 
erogeneous  melting,  the  correlations  in  atomic  positions  are 
retained  on  the  scale  of  the  crystalline  regions  and  the  dif¬ 
fraction  pattern  is  generated  by  superposition  of  the  diffrac¬ 
tion  from  the  liquid  and  crystalline  parts  of  the  system.  In 
other  words,  the  long-range  order  is  more  subtle  in  homoge¬ 
neous  melting  than  in  heterogeneous  melting  due  to  the 
smaller  characteristic  length  scales  at  which  the  homoge¬ 
neous  phase  transformation  takes  place. 

C.  Splitting  of  the  peaks:  The  uniaxial  lattice  expansion 

In  this  section  we  discuss  the  splitting  of  the  diffraction 
peaks  that  takes  place  shortly  after  the  laser  irradiation  [Figs. 
4(b),  5(b),  and  6(b)]  and  demonstrate  that  the  splitting  is  a 
direct  consequence  of  the  uniaxial  thermoelastic  deformation 
of  the  film  in  response  to  the  laser  heating. 

Laser  excitation  of  the  conduction  band  electrons  and  fol¬ 
lowing  electron-phonon  equilibration  lead  to  the  fast  heating 
of  the  lattice  [Fig.  3(a)]  and  generation  of  thermoelastic 
stresses  [Fig.  3(b)],  which  drive  the  expansion  of  the  film. 
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FIG.  7.  (Color  online)  Structure  functions  of  atomic  configura¬ 
tions  generated  by  uniform  uniaxial  deformation  of  a  20  lim  Au  fee 
film  in  the  direction  normal  to  the  free  (001)  surfaces.  The  values  of 
the  deformation  are  indicated  in  the  figure.  Temperature  of  the  film 
is  300  K.  Curves  are  shifted  vertically  with  respect  to  each  other  in 
order  to  better  show  the  changes  in  the  structure  function.  Structure 
peaks  are  identified  by  Miller  indices. 

The  periodic  boundary  conditions  applied  in  the  directions 
parallel  to  the  surfaces  of  the  film  only  allow  the  expansion 
of  the  film  to  proceed  in  the  direction  normal  to  the  surface. 
These  conditions  of  the  lateral  confinement  are  also  realized 
in  experiments,  where  the  laser  spot  diameter  is  much  larger 
then  the  depth  of  the  heated  region  or  thickness  of  the  irra¬ 
diated  film.  The  simulations  are  performed  for  a  single¬ 
crystalline  fee  film  with  (001)  free  surfaces.  The  uniaxial 
deformation  of  the  fee  lattice  along  the  [001]  direction 
changes  the  space  group  symmetry  of  the  lattice  as  the  cubic 
lattice  transforms  into  a  tetragonal  lattice.  Thus,  one  should 
expect  the  occurrence  of  new  diffraction  peaks  correspond¬ 
ing  to  the  face-centered  tetragonal  (fet)  lattice. 

In  order  to  investigate  the  degree  to  which  the  lattice  ex¬ 
pansion  alone  can  explain  the  splittings  observed  in  the  dif¬ 
fraction  profiles  in  Figs.  4(b),  5(b),  and  6(b),  we  calculate 
structure  functions  for  a  series  of  Au  fee  structures  uniaxially 
deformed  along  the  [001]  direction.  The  results  of  the  dif¬ 
fraction  profile  calculations  are  shown  in  Fig.  7.  Splittings 
and  shifts  of  the  diffraction  peaks,  increasing  with  increasing 
deformation,  are  apparent  in  the  figure.  The  diffraction  peaks 
from  (111)  and  (311)  atomic  planes,  present  in  the  original 
fee  structure,  shift  in  the  direction  of  smaller  Q.  New  diffrac¬ 
tion  peaks  with  (002),  (202),  (113),  and  (004)  Miller  indices 
appear  as  the  lattice  transforms  from  fee  to  fct. 

It  should  be  noted  that  the  appearance  of  new  peaks  is  not 
observed  in  MD  simulations  of  a  fee  crystal  slowly  heated  up 
to  the  melting  temperature  under  constant  hydrostatic  pres¬ 
sure  conditions  (with  periodic  boundary  conditions  applied 
in  all  directions).  In  these  simulations  the  effects  of  an  in¬ 
crease  in  temperature  are  limited  to  the  decrease  of  the 
heights  of  the  peaks  due  to  the  atomic  thermal  vibrations  (see 


Sec.  IV  D)  and  shift  of  the  peak  positions  to  smaller  values 
of  Q  due  to  the  isotropic  thermal  expansion  of  the  crystal. 

By  comparing  the  diffraction  profiles  obtained  for  the 
uniaxially  deformed  films  (Fig.  7)  with  the  results  obtained 
in  the  laser  melting  simulations  [Figs.  4(b),  5(b),  and  6(b)], 
we  can  conclude  that  both  the  shifts  and  splittings  of  the 
peaks  can  be  explained  by  the  uniaxial  thermoelastic  expan¬ 
sion  of  the  lattice.  Using  the  (002)  peak  splittings  shown  in 
Fig.  7  as  a  reference,  we  can  estimate  that  uniaxial  deforma¬ 
tions  along  the  (001)  direction  of  2.8%,  5.3%,  and  6.6% 
would  produce  the  values  of  peak  shifts  and  splittings  ob¬ 
served  in  a  simulation  performed  at  180  J/m2  [Fig.  4(b)]  at 
7,  10,  and  1 1  ps,  respectively.  These  values  calculated  from 
the  diffraction  spectra  are  consistent  with  the  ones  obtained 
by  directly  measuring  the  thickness  of  the  film  in  the  snap¬ 
shots  of  the  atomic  configurations  [Fig.  4(a)]  as  well  as  with 
the  distances  between  the  density  peaks  in  the  density  distri¬ 
bution  along  the  direction  normal  to  the  film  surfaces.  The 
distances  between  the  density  peaks  correspond  to  the  spac¬ 
ing  between  the  (001)  lattice  planes  in  the  expanded  film  and 
are  directly  related  to  the  position  of  the  (002)  diffraction 
peak. 

In  the  case  of  the  absorbed  fluence  of  45  J/m2  (Fig.  6), 
the  splittings  of  the  diffraction  peaks  at  10,  60,  180,  and 
500  ps  correspond  to  2.3%,  3.9%,  3.1%,  and  3.1%  uniaxial 
deformations  of  the  lattice.  These  values  are  again  consistent 
with  interplane  distances  measured  in  the  density  distribution 
along  the  direction  normal  to  the  film  surfaces.  The  direct 
measurement  of  deformation  based  on  the  thickness  of  the 
film  is  hampered  in  this  case  by  melting  of  the  surface  re¬ 
gions.  The  diffraction  profiles  calculated  for  times  of  180  and 
500  ps  show  some  subtle  reverse  shifts  and  decrease  in  the 
splittings  of  the  peaks  with  respect  to  the  shift  values  ob¬ 
served  at  60  ps.  These  changes  may  be  related  to  the  gradual 
cooling  of  the  film  associated  with  the  melting  process,  by 
—50  K  from  60  to  180  ps. 

Thus,  the  results  discussed  above  indicate  that,  in  the  flu¬ 
ence  range  considered  in  this  study,  the  melting  of  Au  film  is 
preceded  by  a  significant  thermoelastic  uniaxial  lattice  ex¬ 
pansion,  which  can  be  identified  from  shifts  and  splittings  of 
the  diffraction  peaks.  It  has  been  shown  that  the  uniaxial 
expansion  and  associated  anisotropic  lattice  distortions  can 
significantly  reduce  the  lattice  stability  against  the  initiation 
of  melting  and  can  lead  to  the  homogeneous  nucleation  of 
liquid  regions  at  temperatures  close  to  the  equilibrium  melt¬ 
ing  temperature.20  Although  the  analysis  performed  in  this 
paper  is  done  for  a  single-crystal  film  oriented  perpendicular 
to  [001]  direction,  quantitative  analysis  of  laser-induced  de¬ 
formations  is  also  possible  for  polycrystalline  samples,  as 
soon  as  the  texture  of  the  sample  is  known,  e.g.,  Ref.  48. 
Investigation  of  the  development  of  thermoelastic  deforma¬ 
tions  in  time-resolved  x-ray  or  electron  diffraction  experi¬ 
ments  has  a  potential  for  providing  important  information  on 
the  characteristic  time  scale  of  lattice  heating  and  ther¬ 
moelastic  deformation  as  well  as  on  the  role  of  the  uniaxial 
deformation  in  laser-induced  phase  transformations.  Indeed, 
the  evolution  of  the  lattice  deformation  in  a  Au(lll)  single 
crystal  following  short-pulse  laser  heating  has  been  mea¬ 
sured  with  —10  ps  temporal  resolution  in  x-ray  diffraction 
experiments  and  related  to  the  kinetics  of  the  lattice  tempera- 
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ture  evolution  in  the  surface  region  of  the  irradiated  crystal.49 
Periodic  oscillations  of  the  diffraction  peak  positions  have 
been  recently  probed  with  —0.5  ps  resolution  in  electron  dif¬ 
fraction  experiments  and  related  to  the  elastic  vibrations  of  a 
free  standing  A1  film  irradiated  with  a  femtosecond  laser 
pulse.15,16  Similar  oscillations  of  the  diffraction  peak  posi¬ 
tions  have  been  observed  in  simulations  of  200  fs  pulse  laser 
excitation  of  a  Ni  film  at  low  fluences,  below  the  threshold 
for  laser  melting.50  A  discussion  of  the  implications  of  the 
computational  predictions  on  the  shifts  and  splittings  of  the 
diffraction  peaks  for  interpretation  of  experimental  results 
obtained  for  polycrystalline  targets  and  a  fixed  angle  of  inci¬ 
dence  of  the  probe  beam  to  the  target  surface  is  given  in  Sec. 
V. 


D.  Thermal  effects 


The  decrease  of  the  intensity  of  the  diffraction  peaks  in 
the  irradiated  films  can  result  from  both  structural  changes 
and  increasing  thermal  vibrations  of  the  atoms.  The  latter 
contribution,  related  to  the  thermal  smearing-out  of  the  lat¬ 
tice  planes,  can  be  described  by  the  Debye-Waller  factor 
e~2M,  which  relates  the  decrease  of  the  peak  intensity  to  the 
temperature  evolution  in  the  system,35 


S(G,')-I=[s(e,0ps)-1] 


exp[-  2 M(T(t),Q)] 
exp[-  2M(T(0  ps),Q)] 


(ID 


The  temperature  dependence  of  the  quantity  M  in  the  Debye- 
Waller  factor  can  be  expressed  through  the  mean-square  dis¬ 
placement  (MSD)  (u2)  of  an  atom  vibrating  around  its  equi¬ 
librium  lattice  position, 


e-2M  =  e 


^^sin2  0(m2)  (m2)Q2 


3\ 


=  e 


(12) 


The  MSD  can  be  computed  directly  from  atomic  trajecto¬ 
ries  obtained  in  constant-volume  MD  simulations,  performed 
for  a  range  of  temperatures  corresponding  to  those  measured 
in  the  simulations  for  laser  excitation.  The  results  of  such 
calculations  can  be  then  used  to  predict  the  reduction  of  the 
peak  heights  in  the  laser  excitation  simulations  due  to  the 
increase  of  the  lattice  temperature  alone. 

The  results  of  the  application  of  this  approach  to  the  simu¬ 
lation  performed  at  an  absorbed  fluence  of  180  J/m2  are 
shown  in  Fig.  8  for  (111),  (220),  and  (311)  diffraction  peaks. 
To  make  the  connection  between  the  thermal  effects  de¬ 
scribed  by  Eqs.  (11)  and  (12)  and  the  laser  excitation  simu¬ 
lations,  we  first  calculate  the  average  temperatures  of  the  film 
at  different  times.  The  temperatures  are  then  used  to  calcu¬ 
late  the  MSD  and  corresponding  reduction  of  the  peaks  from 
the  Debye- Waller  factor,  Eqs.  (11)  and  (12).  The  results  of 
the  calculations  are  shown  in  Fig.  8  by  dashed  lines. 

The  heights  of  the  peaks  in  the  diffraction  spectra  calcu¬ 
lated  for  the  laser  excitation  simulation  are  shown  in  Fig.  8 
up  to  14  ps,whereas  the  results  of  the  Debye- Waller  calcula¬ 
tions  are  shown  up  to  9  ps,  the  time  when  the  average 
temperature  of  the  him  reaches  1 ,25Tm  =  1203  K  in  the  simu¬ 
lation  of  laser  melting  (above  this  temperature  a  fast  homo¬ 
geneous  melting  takes  place  within  several  picoseconds 19-21 
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FIG.  8.  (Color  online)  Normalized  height  of  the  (111),  (220), 
and  (311)  peaks  as  a  function  of  time  after  irradiation  of  a  20  nm 
Au  Him  with  a  200  fs  laser  pulse  at  an  absorbed  fluence  of 
180  J/m2.  All  peak  heights  are  normalized  to  their  values  at  0  ps 
(T=300  K).  Dotted  lines  are  calculated  through  the  Debye- Waller 
factor  up  to  the  temperature  of  1.25  Tm,  which  corresponds  to  the 
limit  of  crystal  stability.  The  values  of  the  average  temperature  of 
the  film  at  different  times  after  the  laser  pulse  are  shown  in  the 
additional  top  x  axis. 

and  the  MSD  for  atomic  vibration  in  a  crystal  cannot  be 
obtained).  There  is  a  good  agreement  between  the  normal¬ 
ized  peak  heights  and  the  Debye- Waller  calculations  up  to 
—4  ps,when  the  average  lattice  temperature  is  —750  K.  This 
agreement  indicates  that  the  observed  reduction  of  the  dif¬ 
fraction  peaks  during  the  first  4  ps  is  primarily  due  to  the 
increased  thermal  vibration  of  the  atoms.  Starting  from 
5  ps, significant  deviations  from  the  Debye- Waller  calcula¬ 
tions  can  be  observed  for  (220)  and  (311)  peaks.  As  shown  in 
Figs.  2,  3,  and  4(a),  the  melting  process  starts  only  at  — 10  ps 
in  this  simulation  and  the  onset  of  the  deviation  of  the  peak 
intensities  from  the  Debye-Waller  calculations  can  be  attrib¬ 
uted  to  the  uniaxial  film  expansion  discussed  in  Sec.  IV  C. 
As  we  can  see  from  Fig.  3(b), the  accumulated  compressive 
pressure  increases  during  the  first  picoseconds  after  the  laser 
pulse  and  drives  the  expansion  of  the  film.  The  expansion 
results  in  the  distortion  of  the  fee  lattice  and  associated  te¬ 
tragonal  splitting  of  the  (220)  from  the  (202)/(022)  and  the 
(113)  from  the  (311)/(131)  peaks.  The  splitting  is  apparent  in 
the  structure  function  shown  in  Fig.  4(b)  for  7  ps,  and  it 
starts  to  contribute  to  the  decrease  of  the  intensity  of  indi¬ 
vidual  peaks  from  —4.5  ps.  The  (111)  peak  does  not  split  and 
the  data  points  are  well  described  by  the  Debye- Waller  cal¬ 
culations  up  to  the  onset  of  melting  at  about  9  ps  (Fig.  8). 
During  the  melting  process  the  (111)  peak  starts  to  overlap 
with  a  broad  first  peak  characteristic  of  the  liquid  structure 
and  the  points  for  the  (111)  peak  plotted  in  Fig.  8  for  times 
starting  from  10  ps  correspond  to  the  heights  measured  from 
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the  background  level  provided  by  the  broad  liquid  peak  of 
the  structure  function  (Fig.  4).  A  sharp  drop  of  the  intensity 
of  the  (111)  peak  is  observed  during  the  melting  process, 
from  10  to  14  ps. 

Thus,  we  can  conclude  that  the  reduction  of  the  diffrac¬ 
tion  peaks  heights  is  mainly  defined  by  the  increasing  lattice 
temperature  (Debye- Waller  factor)  during  the  first  several  pi¬ 
coseconds  following  the  laser  excitation,  when  the  film  does 
not  have  time  to  expand  and  is  “inertially  confined.”22  As  the 
film  expands  in  response  to  the  laser-induced  thermoelastic 
stresses,  the  peak  splitting  due  to  the  uniaxial  lattice  defor¬ 
mations  starts  to  contribute  to  the  reduction  of  the  peaks. 
Finally,  during  the  melting  process  the  reduction  of  the 
height  of  the  diffraction  peaks  is  dominated  by  the  destruc¬ 
tion  of  the  crystal  order. 

A  similar  analysis  of  the  effect  of  the  increasing  lattice 
temperature  on  the  diffraction  peak  intensities  measured  in 
the  time-resolved  electron  diffraction  study  of  laser-driven 
melting  of  thin  A1  films  is  presented  in  Ref.  12.  The  Debye 
equation  and  TTM  are  used  in  this  work  to  predict  the  atomic 
mean  square  displacements  and  the  temperature  evolution, 
respectively.  The  conclusions  from  the  analysis  of  the  experi¬ 
mental  data  are  in  a  good  qualitative  agreement  with  the 
computational  results  discussed  above.  The  reduction  of  the 
diffraction  peak  heights  initially  follows  the  Debye- Waller 
calculation,  with  the  onset  of  the  deviations  attributed  to  the 
loss  of  the  crystalline  order.  A  much  shorter,  as  compared  to 
the  simulations,  time  for  the  beginning  of  the  melting  pro¬ 
cess,  1.5  ps,  can  be  explained  by  a  much  stronger  electron- 
phonon  coupling  in  A1  as  compared  to  Au.  Indeed,  in  a  simu¬ 
lation  performed  for  an  A1  film  at  irradiation  conditions 
comparable  to  the  experimental  ones,  we  obtain  an  excellent 
quantitative  agreement  in  both  the  time  of  the  onset  of  melt¬ 
ing  and  the  duration  of  the  melting  process  (see  Sec.  V). 

E.  Real-space  correlations 

Structural  changes  during  the  laser-induced  phase  trans¬ 
formations  can  be  further  investigated  by  keeping  track  of 
the  evolution  of  real-space  correlation  functions  calculated 
for  atomic  configurations  predicted  in  MD  simulations.  In 
particular,  correlations  in  atomic  positions  can  be  expressed 
in  the  form  of  the  reduced  pair  distribution  function  G(r), 
which  describes  the  deviation  of  the  pair  density  function 
p{r),  defined  by  Eq.  (6),  from  the  average  density  p0, 

G(r)  =  4Trr[p(r)  -  p0\.  (13) 

Experimentally,  G(r)  can  be  calculated  by  the  sine  Fourier 
transform  of  the  structure  function  S(Q)  and  requires  the 
ability  to  obtain  the  diffraction  intensity  of  an  adequate  qual¬ 
ity  over  a  broad  range  of  Q.  This  requirement  presents  a 
challenge  for  ultrafast  time-resolved  diffraction  experiments 
that  are  often  limited  to  the  analysis  of  the  diffraction  signals 
obtained  for  fixed  Bragg  angles.  The  challenge  of  performing 
measurements  for  a  range  of  scattering  vectors  can  be  met  by 
increasing  intensity  of  the  x-ray  or  electron  probes  and  accu¬ 
mulating  the  diffraction  intensity  over  a  large  number  of 
pulses. 11,12’14,15  In  simulations,  both  p(r)  and  G(r)  can  be 
calculated  directly  from  atomic  configurations,  using  Eqs.  (6) 
and  (13). 


An  important  advantage  of  G{r )  is  that  its  correlation 
peaks  provide  direct  information  on  the  structural  coherence 
in  the  system.  In  real  crystals,  G(r)  obtained  from  diffraction 
experiments  oscillate  around  zero  with  peak  amplitudes 
gradually  decreasing  with  increasing  r  due  to  the  crystal  im¬ 
perfections  and  the  finite  resolution  of  Q  value  in  the 
measurement.38  For  disordered  (liquid  or  amorphous)  struc¬ 
tures,  the  amplitude  of  G(r)  oscillations  decreases  much 
faster,  indicating  the  absence  of  long-range  order.51 

In  Fig.  9  the  reduced  pair  distribution  function  G(r)  is 
presented  for  a  simulation  of  a  20  nm  Au  film  irradiated  with 
a  200  fs  laser  pulse  at  an  absorbed  fluence  of  180  J/m2.  Be¬ 
fore  the  laser  excitation  (0  ps),  G(r)  exhibits  the  characteris¬ 
tic  features  of  the  fee  crystalline  structure.  Each  peak  in  G(r) 
corresponds  to  a  specific  interatomic  distance  between  a  pair 
of  atoms  in  a  perfect  fee  structure.  For  example,  the  first 
peak  in  G(r)  matches  the  distance  between  the  nearest  neigh¬ 
bors  (2.89  A)  in  a  fee  crystal  equilibrated  at  300  K.  Changes 
of  G(r)  during  the  fast  homogeneous  melting  (see  Fig.  4)  are 
shown  in  Fig.  9(b).  By  comparing  G(r )  at  0  and  1 1  ps  we 
can  see  that  before  the  onset  of  melting  the  increasing  lattice 
temperature  results  in  a  significant  broadening  of  the  corre¬ 
lation  peaks  and  reduction  of  their  intensities.  Some  of  the 
peaks  become  completely  obscured  and  merge  with  the 
neighboring  ones,  e.g.,  the  peak  corresponding  to  the  second- 
neighbor  shell  cannot  be  identified  in  Fig.  9(b).  The  seeming 
“disappearance”  of  the  correlation  peaks  before  the  onset  of 
melting,  however,  does  not  correspond  to  any  structural 
changes  in  the  crystal  but  is  just  a  consequence  of  the  ther¬ 
mal  broadening  of  all  peaks,  as  demonstrated  in  the  analysis 
presented  at  the  end  of  this  section.  Despite  the  broadening, 
the  peaks  at  11  ps  remain  well  defined  for  the  range  of  in¬ 
teratomic  distances  considered  in  the  calculation,  confirming 
the  presence  of  the  long-range  crystalline  ordering  in  the 
atomic  configuration.  Between  the  times  of  11  and  14  ps  the 
peaks  in  the  high-r  region  disappear,  reflecting  a  complete 
loss  of  the  long-range  order  in  the  solid-to-liquid  transition. 
The  remained  correlation  peaks  continue  to  broaden  and  di¬ 
minish  between  the  times  of  14  and  40  ps  [Fig.  9(c)],  indi¬ 
cating  that  there  are  still  some  substantial  changes  in  the 
short-  or  medium-range  correlations  in  atomic  positions.  This 
observation  suggests  that  within  this  period  of  time  the  liquid 
structure  created  after  the  collapse  of  crystalline  state  is  still 
in  a  nonequilibrium  state  and  retain  some  “memory”  of  the 
crystalline  state.52  The  continuing  increase  of  the  lattice  tem¬ 
perature  [Fig.  3(a)]  facilitates  the  destruction  of  the  remain¬ 
ing  medium-range  ordering. 

An  alternative  way  of  representing  real-space  correlations 
in  atomic  positions  is  provided  by  the  radial  distribution 
function  (RDF)  R(r)  =  4Trr2p(r),  from  which  coordination 
numbers  corresponding  to  different  neighbor  shells  can  be 
determined  by  calculating  the  surface  areas  under  the  corre¬ 
sponding  correlation  peaks.  The  surface  area  can  be  found  by 
the  integration  of  R(r)  with  the  integration  limits  correspond¬ 
ing  to  the  local  minima  on  the  sides  of  the  corresponding 
peak  of  the  RDF.  The  RDF  obtained  in  a  simulation  per¬ 
formed  for  a  20  nm  Au  film  irradiated  by  a  200  fs  laser  pulse 
at  180  J/m2  is  shown  for  0,  10,  and  20  ps  after  the  laser 
pulse  in  Fig.  10.  For  0  ps,  the  first  and  second  coordination 
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FIG.  9.  (Color  online)  Reduced  pair  distribution  function  G(r)  computed  for  atomic  configurations  obtained  in  a  simulation  of  a  20  nm 
Au  film  irradiated  with  a  200  fs  laser  pulse  at  an  absorbed  fluence  of  180  J/m2  at  times  of  (a)  0,  (b)  11-14,  and  (c)  14  and  40  ps. 


numbers  corresponding  to  the  first  and  second  peaks  of  the 
RDF  are  calculated  to  be  11.92  and  5.96.  The  values  are  less 
than  the  coordination  numbers  in  a  perfect  fee  lattice,  12  and 
6,  due  to  the  presence  of  the  two  free  surfaces  in  our  system, 
where  atoms  have  a  smaller  number  of  neighbors.  At  10  ps 
the  first  coordination  number  is  determined  to  be  12.1  ±0.4 
with  the  error  value  related  to  the  uncertainty  in  defining  the 
local  minima  in  Fig.  10.  This  result  indicates  that,  in  agree¬ 
ment  with  visual  picture  of  atomic  configuration  at  10  ps  in 
Fig.  4(a),  the  local  crystalline  structure  at  this  time  is  still 
preserved  and  most  atoms  still  maintain  12  nearest  neigh¬ 
bors.  By  the  time  of  20  ps  the  first  coordination  number  de¬ 
creases  down  to  10.9 ±0.4,  reflecting  the  destruction  of  the 
crystalline  structure  that  is  completed  by  this  time  [Fig.  4(a)]. 
A  similar  decrease  of  the  first  coordination  number  from  12 
to  10  during  the  laser-induced  melting  process  has  been  de¬ 
duced  from  time-resolved  electron  diffraction  measurements 
reported  in  Ref.  1 1  for  a  20  nm  A1  film. 


The  output  of  MD  simulations  contains  complete  infor¬ 
mation  on  the  atomic  dynamics  induced  by  laser  excitation 
and  allows  for  detailed  investigation  of  the  origin  of  the 
changes  in  the  correlation  functions  discussed  above  and  il¬ 
lustrated  in  Figs.  9  and  10.  In  particular,  one  question  arises 
as  to  why  the  correlation  peak  corresponding  to  the  second- 
neighbor  shell  in  the  fee  structure  [Fig.  9(a)]  completely  dis¬ 
appears  even  before  the  onset  of  the  melting  process,  e.g.,  the 
plots  for  11  and  10  ps  in  Figs.  9  and  10,  respectively.  To 
answer  this  question  and  to  better  understand  other  aspects  of 
the  evolution  of  the  correlation  functions,  we  trace  the  mo¬ 
tions  of  the  first-,  second-,  and  third-nearest  neighbors  of 
each  atom  in  the  original  fee  lattice.  By  doing  so,  we  decom¬ 
pose  the  RDF  into  separate  contributions  from  the  first  three 
original  coordination  shells.  The  result  of  the  decomposition, 
shown  in  Fig.  11,  indicates  that  before  10  ps,  on  average, 
most  atoms  in  the  three  coordination  shells  are  still  localized 
around  their  equilibrium  positions.  Significant  broadening  of 
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Distance  r  (A) 

FIG.  10.  (Color  online)  Radial  distribution  function  R(r)  com¬ 
puted  for  atomic  configurations  obtained  in  a  simulation  of  a  20  nm 
Au  film  irradiated  with  a  200  fs  laser  pulse  at  an  absorbed  fluence 
of  180  J/m2  at  times  of  0,  10,  and  20  ps.  Coordination  number  Nc 
for  the  first  coordination  shell  at  ~2.89  A  can  be  obtained  by  inte¬ 
gration  of  R(r)  between  the  two  local  minima  on  either  side  of  the 
first  peak. 

all  the  peaks,  however,  completely  obscures  the  contribution 
from  the  peak  corresponding  to  the  second-neighbor  shell 
[Fig.  11(d)],  which  shows  itself  only  as  a  shoulder  on  the 
new  second  peak  located  at  the  position  corresponding  to  the 
third-neighbor  shell. 

After  the  collapse  of  crystalline  order,  atoms  gain  the  free¬ 
dom  to  diffuse  around.  For  the  first-nearest-neighbor  shell, 
atoms  can  only  move  outward,  which  leads  to  the  long  tail  in 
the  original  first-nearest-shell  contribution  to  the  RDF  at  14 
and  20  ps  [Fig.  11(a)].  The  gradual  disappearance  of  the 
original  peak  can  be  explained  by  the  high  mobility  of  atoms 
in  the  liquid  state,  leading  to  the  loss  of  correlation  between 
the  positions  of  the  nearest  neighbors  in  the  original  fee  lat¬ 
tice.  A  similar  gradual  disappearance  of  the  peak  correspond¬ 
ing  to  the  original  third-nearest-neighbor  shell  is  observed  in 
Fig.  11(c),  with  some  of  the  neighbors  moving  closer  to  the 
atoms  of  origin,  others  moving  further  away  from  them. 

The  evolution  of  the  peak  corresponding  to  the  original 
second-nearest-neighbor  shell  is  rather  different  from  the 
ones  of  the  first  and  third  peaks  [Fig.  11(b)].  By  the  time  of 
14  ps  the  second  fee  peak  completely  disappears,  splitting  its 
intensity  between  the  two  peaks  present  in  the  liquid  state. 
The  fast  disappearance  of  the  second  peak,  also  observed  in 
a  recent  time-resolved  electron  diffraction  study,11  suggests 
that  the  interatomic  distance  corresponding  to  the  second- 
nearest-neighbor  shell  in  the  fee  structure  is  not  characteristic 
of  the  short-range  order  in  the  liquid  structure.  This  inter¬ 
atomic  distance  appears  in  the  octahedral  atomic  configura¬ 
tion  characteristic  of  the  close-packed  crystals  but  is  not  typi¬ 
cal  for  tetrahedral  clusters  responsible  for  the  short-range 
order  in  one-component  noncrystalline  structures.53  Thus,  as 
soon  as  material  melts,  the  atoms  tend  to  escape  from  this 


interatomic  distance — octahedral  configurations  disappear 
and  a  broad  second  peak  corresponding  to  the  interatomic 
distances  found  in  clusters  with  local  tetrahedral  and  icosa- 
hedral  short-range  order54,55  develops. 

V.  CONCLUSION  AND  CONNECTIONS  TO  EXPERIMENTS 

Computational  analysis  of  the  structure  function  and  pair 
distribution  function,  performed  for  transient  atomic  configu¬ 
rations  generated  in  MD  simulations  of  fast  laser-induced 
phase  transformations,  provides  a  connection  between  the 
atomic-level  structural  rearrangements  and  changes  in  the 
diffraction  profiles.  The  laser-induced  evolution  of  the  dif¬ 
fraction  peaks  is  found  to  be  defined  by  the  combination  of 
the  following  three  factors:  (1)  increasing  amplitude  of 
atomic  vibrations  associated  with  the  fast  temperature  rise, 
(2)  peak  splitting  due  to  the  uniaxial  lattice  deformations, 
and  (3)  destruction  of  the  crystal  order  in  the  melting  pro¬ 
cess.  The  contribution  of  the  increasing  thermal  atomic  vi¬ 
brations  dominates  during  the  first  several  picoseconds  after 
the  laser  pulse  and  can  be  well  described  by  the  Debye- 
Waller  factor.  The  film  expansion  in  response  to  the  laser- 
induced  thermoelastic  stresses  results  in  shifts  and  splittings 
of  the  diffraction  peaks,  providing  an  opportunity  for  experi¬ 
mental  probing  of  the  ultrafast  deformations.  Finally,  the  on¬ 
set  of  the  melting  process  results  in  further  reduction  the 
diffraction  peaks  due  to  the  destruction  of  the  crystal  order. 

Two  of  the  three  contributions  to  the  reduction  of  the 
diffraction  peak  intensity  discussed  above  have  been  ob¬ 
served  in  recent  time-resolved  electron  diffraction 
experiments. 11-13  For  a  20  nm  A1  film  irradiated  at  an  exci¬ 
tation  fluence  of  700  J/m2,  the  Debye- Waller  calculations 
are  found  to  be  in  a  good  agreement  with  experimental  data 
for  the  first  1.5  ps,  whereas  further  reduction  of  the  height  of 
the  diffraction  peaks  has  been  attributed  to  the  ultrafast  melt¬ 
ing  that  completes  by  the  time  of  3.5  ps  after  the  laser  pulse. 
No  shifts  in  the  positions  of  the  diffraction  peaks  and  no 
splitting  of  the  peaks  have  been  observed  in  this  study  as 
well  as  in  time-resolved  x-ray  diffraction  investigations  of 
even  faster  nonthermal  melting  processes.5,6  The  absence  of 
the  shifts  and  splittings  of  the  diffraction  peaks  can  be  related 
to  the  condition  of  the  inertial  stress  confinement,22,56  when 
the  lattice  heating  and  melting  are  taking  place  as  constant- 
volume  processes.  In  the  case  of  the  thermal  melting,  the 
time  of  the  lattice  heating  is  defined  by  the  laser  pulse  dura¬ 
tion  Tp,  and  the  time  of  the  electron-phonon  equilibration, 
re- ph,  whichever  is  larger.  The  condition  for  the  inertial  stress 
confinement  can  be  then  formulated  as  max{  rp ,  Te.ph}  =S  r, 
~d/Cs,  where  d  is  the  film  thickness  and  Cs  is  the  speed  of 
sound  in  the  film  material.  Due  to  the  strong  electron-phonon 
coupling  in  Al,  the  time  of  the  electron-phonon  equilibration 
(defined  by  the  exponential  fit  of  the  time  dependence  of  the 
energy  transferred  from  the  electrons  to  the  lattice)  is  short, 
Te_ ph~  1.5  ps,  and  fast  homogeneous  melting  takes  place  be¬ 
fore  any  significant  expansion  of  the  film  can  take  place. 
Indeed,  the  results  of  the  simulations  of  laser  melting  of  a 
20  nm  Al  film  performed  for  experimental  conditions  re¬ 
ported  in  Ref.  11  show  an  ultrafast  melting  process  within 
the  first  3-4  ps  with  no  peak  shifts  and  splittings  observed  in 
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FIG.  11.  (Color  online)  Decomposed  radial  distribution  function  R(r)  of  (a)  first,  (b)  second,  and  (c)  third  coordination  shells  computed 
for  atomic  configurations  obtained  in  a  simulation  of  a  20  nm  Au  film  irradiated  with  a  200  fs  laser  pulse  at  an  absorbed  fluence  of  180  J/m2 
at  times  of  0,  10,  14,  and  20  ps.  The  peaks  are  calculated  for  atoms  that  belong  to  the  corresponding  coordination  shells  at  the  time  of  0  ps. 
The  data  for  the  three  coordination  shells  and  their  sum  at  10  ps  are  shown  in  (d). 


the  diffraction  spectra  (Fig.  12).  The  time  of  the  melting 
onset,  1 .5  ps,  the  duration  of  melting  process,  and  the  evolu¬ 
tion  of  the  structure  function  are  all  in  a  very  good  quantita¬ 
tive  agreement  with  experimental  observations.  In  the  case  of 
even  more  rapid  nonthermal  melting  occurring  on  a  subpico¬ 
second  time  scale5-6  there  is  no  doubt  that  the  destruction  of 
the  crystal  order  occurs  before  any  expansion  of  the  material 
can  take  place. 

In  the  simulations  reported  in  this  paper  for  Au  films,  a 
weak  electron-phonon  coupling  in  Au  leads  to  a  relatively 
slow  heating  of  the  film,  re_ph~  15  ps,  whereas  the  small 
thickness  of  the  film  allows  for  a  fast  relaxation  of  ther¬ 
moelastic  stresses  with  ts~  10  ps.  Actually,  due  to  the  pres¬ 
ence  of  two  free  surfaces,  the  relaxation  of  the  initial  com¬ 
pressive  pressure  in  a  freestanding  film  takes  place  at  an  even 


shorter  time  scale  of  ~5  ps  and  is  defined  by  the  propagation 
of  the  two  unloading  waves  from  the  free  surfaces  of  the  film 
[Fig.  3(b)].  As  a  result,  the  uniaxial  expansion  of  the  film 
precedes  the  onset  of  the  structural  transformation,  and  shifts 
and  splittings  of  the  diffraction  peaks  are  observed  in  all 
simulations  performed  for  Au  films. 

A  clear  separation  of  the  time  scales  for  lattice  heating 
and  melting  has  also  been  observed  in  the  first  time-resolved 
electron  diffraction  experiments  performed  for  20  nm  Au 
films.13  At  an  absorbed  laser  fluence  of  1 19  J/m2,  well  above 
the  threshold  for  the  complete  melting  of  the  film,  the  melt¬ 
ing  starts  at  about  7  ps  after  the  laser  pulse  and  is  completed 
in  about  3  ps.  While  the  duration  of  the  melting  process  is  in 
an  excellent  agreement  with  the  results  of  our  simulations 
(Fig.  2),  the  time  of  the  melting  onset  is  significantly  shorter 
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FIG.  12.  (Color  online)  Structure  functions  calculated  for  atomic  configurations  obtained  in  a  simulation  of  a  20  nm  A1  film  irradiated 
with  a  120  fs  laser  pulse  at  an  absorbed  fluence  of  84  J/m2  at  times  of  (a)  0,  (b)  1,  (c)  2,  and  (d)  3  ps.  Insets  show  the  corresponding 
snapshots  of  the  atomic  configurations,  with  the  direction  of  the  laser  pulse  shown  in  (a).  Atoms  are  colored  according  to  the  local  order 
parameter  (blue  atoms  have  local  crystalline  surroundings;  red  atoms  belong  to  the  liquid  phase).  The  irradiation  conditions  are  the  same  as 
in  the  experiment  reported  in  Ref.  11  (incident  fluence  of  700  J/m2,  reflectivity  of  88%). 


in  the  experiment  as  compared  to  the  simulations,  7  ps  in 
experiment  vs  17  ps  in  a  simulation  performed  for  compa¬ 
rable  irradiation  conditions.  This  quantitative  discrepancy 
can  be  attributed  to  the  assumption  of  the  temperature- 
independent  electron-phonon  coupling  constant  used  in  this 
work.  Indeed,  recent  calculations  accounting  for  the  contri¬ 
bution  of  the  thermal  excitation  of  the  d  band  electrons  to  the 
electron  heat  capacity  and  electron-phonon  coupling57  sig¬ 
nificantly  increase  the  rate  of  the  lattice  heating  and  show  an 
excellent  agreement  with  experimental  results. 

Note  that  the  structure  functions  discussed  in  this  paper 
(Figs.  4-7  and  12)  are  for  single  crystals  and  include  the 


information  on  all  the  lattice  plane  systems  present  in  the 
crystal.  In  polycrystalline  samples  different  orientations  of 
crystallographic  planes  will  be  affected  in  a  different  way  by 
the  uniaxial  expansion  of  the  film  and  the  structure  function 
would  exhibit  broadening  of  the  peaks  instead  of  the  split¬ 
ting.  In  experiments  performed  with  a  fixed  angle  of  inci¬ 
dence  of  the  probe  beam  to  the  target  surface,  only  the  lattice 
planes  that  have  the  correct  Bragg  angle  with  the  probe  beam 
contribute  to  the  measured  diffraction  spectrum.  Splitting  of 
the  peaks  can  only  be  observed  experimentally  if  probing  is 
done  at  different  incidence  angles,  so  that  different  orienta¬ 
tions  of  the  same  system  of  lattice  planes  with  respect  to  the 
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direction  of  the  uniaxial  thermoelastic  lattice  expansion  is 
probed.  The  maximum  value  of  the  peak  shift  is  defined  in 
this  case  not  only  by  the  extent  of  the  film  expansion  but  also 
by  the  orientation  of  the  corresponding  lattice  planes  with 
respect  to  the  direction  of  the  film  expansion.  The  peak  shift 
can  be  rather  small  if  the  planes  having  a  small  Bragg  angle 
are  probed  in  the  transmission  mode  by  an  electron  beam 
normal  to  the  surface  of  the  film.1316  Nevertheless,  for  a 
given  direction  of  the  probe  beam,  analysis  of  the  peak  shifts 
in  time-resolved  x-ray  or  electron  diffraction  experiments 
can  provide  valuable  information  on  the  characteristic  time 
scale  of  the  lattice  heating  and  thermoelastic 
deformation1516'49  as  well  as  on  the  role  of  the  uniaxial  de¬ 
formation  in  laser-induced  phase  transformations.20  A  discus¬ 
sion  of  the  simulation  results  for  Ni  films,  where  periodic 
oscillations  of  the  diffraction  peak  positions  are  observed 
under  low-fluence  irradiation  conditions  (below  the  threshold 
for  melting)  and  related  to  the  laser-induced  elastic  vibra¬ 
tions  of  the  film,  is  presented  elsewhere.50 

The  emerging  time-resolved  electron  and  x-ray  diffraction 
probe  experiments  are  opening  new  opportunities  for  inves¬ 
tigation  of  atomic  dynamics  in  the  time  and  spatial  domains 
accessible  for  direct  atomistic  simulations.  The  results  of  the 


simulations  provide  a  direct  link  between  the  experimental 
observations  and  atomic-level  structural  changes  in  the  irra¬ 
diated  material  and  help  in  interpretation  of  experimental 
data.  In  particular,  the  observed  differences  in  the  evolution 
of  the  diffraction  profiles  in  the  homogeneous  and  heteroge¬ 
neous  melting  processes  can,  along  with  kinetics  arguments, 
be  instrumental  in  distinguishing  between  the  two  melting 
mechanisms.  The  shift  and  splitting  of  the  diffraction  peaks 
reflect  the  uniaxial  thermoelastic  lattice  deformation  of  the 
film  prior  to  melting  and  can  be  used  for  experimental  prob¬ 
ing  of  the  laser-induced  ultrafast  deformations. 
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